Self-organization of dissipationless solitons in negative refractive index materials 
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General nonlinear and nonparaxial dissipative complex Helmholtz equations for magnetic and 
electric fields propagating in negative refractive index materials (NIMs) are derived ab initio from 
Maxwell equations. In order to describe nonconservative soliton dynamics in NIMs, such coupled 
equations are reduced into generalized Ginzburg-Landau equation. Cross-compensation between the 
excess of saturating nonlinearity, losses, and gain renders these self-organized solitons dissipationless 
and exceptionally robust. The presence of such solitons makes NIMs effectively dissipationless. 



Metamaterials are novel artificial composite structures 
manufactured in order to display various peculiar very 
promising properties [l| . The most studied nowadays are 
negative-refractive-index materials (NIMs) assembled in 
such a way to exhibit simultaneously negative effective 
permittivity e and permeability A number of 

challenging optical NIM devices is supposed to work in 
interaction with lasers As a consequence, nonlinear 
effects during the propagation of electromagnetic (EM) 
waves in NIMs have to be taken into account. In par- 
ticular, in such materials conservative localized solitonic 
structures of Schrodinger type are considered [E]i0i0 
Q. However, in real media dissipation is always present. 
Whenever the absorption is weak, it can be ignored pro- 
vided that EM field frequency lays beyond the resonance 
bands of medium. However, in NIMs the problem of EM 
field absorption is the central issue due to the resonant 
character of interactions. Persistence of losses prevents 
also the generation of Schrodinger type solitons in or- 
dinary materials, a fortiori in NIMs. Numerous highly 
desirable properties of NIMs like the superresolution [Ij] 
and the cloaking [Sjl, are either altered or annihilated 
by the dissipation [9|. A wide class of dissipative sys- 
tems, ranging from nonlinear optics, plasma physics, and 
fluid dynamics to superfluidity, superconductivity, and 
Bose-Einstein condensates, can be modeled by complex 
Ginzburg-Landau equation . 

In order to study propagation of EM field in optical 
NIMs, we derive, in this letter, ab initio from Maxwell 
equations (ME) general nonlinear and nonparaxial dis- 
sipative complex Helmholtz equations for magnetic and 
electric fields. A very promising class of self-organized 
localized EM structures are spatiotemporal solitons ren- 
dered dissipationless due to the cross-compensation of 
the excess of saturating nonlinearity, losses, and gain 11 1. 
Obtained coupled equations are reduced into generalized 
Ginzburg-Landau equation in order to describe for the 
first time dissipationless soliton dynamics in NIMs. Such 
a NIM in presence of dissipationless solitons may be con- 
sidered as an effectively dissipationless novel active com- 
posite metamaterial. 

The propagation of EM radiation in a nonlinear me- 
dia is described by ME rotE = -dB/dt, rotH =dT>/dt 
and constitutive relations between the magnetic induc- 



tion B and the magnetic field H as well as between 
the electric field E and the electric induction D. The 
response of the medium to a quasi-monochromatic EM 
wave is considered. Real vectorial fields F = (E, B, H, D) 
read F(r,t) = F exp (ik^r — wi) -I- c.c, where F is 
slowly varying function in space and time and uj is 
the carrying pulse frequency. The complex wave vec- 
tor h = \ir + iki is determined by the linear dispersion 
relation k^c^/w^ — e {lu) jj. [uj) . Respectively e and /i 
are complex permittivity and permeability of a dissipa- 
tive medium. In NIMs the refractive index is nega- 
tive n = — (3fi[e/Lt])^/^ = krc/oj, hence, the real part of 
wave vector kr is negative too. In the case of tempo- 
ral dispersive media e (w) and ii (w) are expanded in the 
series around the carrying frequency up to the second 
order 0,[lll dTi/dt « e^{a,d^/dt + ia'j2d^E/dt^ - 
iujeE - icje„iE)e~''^* and dB/dt w fil{a^,dli/dt + 
ia'^/2d'^U/dt'^ - iojfiH. ~ itj/i„;H)e-*"*, where = 
d{Loe)/dLu, a'^ = da^jdio, — d{LLi^i)/du}, and a'^ = 
da^/dbj. In order to take into account the nonlinear loss 
and gain the nonlinear permittivity £„( (|Ep) — 5R[e„;]-f 
i9[e„i] and permeability (|Hp) = 3fJ[/i„i] -f 
are considered as complex. Taking curl of ME and ne- 
glecting vectorial terms (V (V • E) — 0) yield following 
equations rewritten in new variables C, = z and t = 



t — z jvg with real group velocity Vg 
(5R[w (/lae + ea^) (2A:c)-i])"\ 



(9 (nuj) jduj) 



and 



2ifcf -f f;i + AE-t-^M£„iE 
-ujHolJ-ni^ X H -f i^^[e^l]I] = 0, 



-w£o£nik X E + i^^[en]li = 0, 



(1) 



(2) 



with the complex function W — cVg' — O.buj^a'^ — 
Q.^ujea'^ — a^j^a^. The second order z derivatives in AE 
and AH should be kept in Eqs. (1-2) not only to ac- 
count for nonparaxial pulse dynamics but also to have an 
adequate description of the strong dissipation. In best 
of our knowledge it is for the first time that such gen- 
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eralized coupled complex Helmholtz equations describ- 
ing simultaneously nonparaxial and dissipative effects in 
NIMs are ab initio derived from ME. However, investiga- 
tion of the EM filed dynamics based on the generalized 
Helmholtz equations is beyond of intended scope of the 
current letter and will be presented elsewhere. There- 
fore, in what follows only weak dissipation is consid- 
ered. As a consequence, the imaginary part of linear 
permittivity e^, permeability /Lti, and wave vector ki are 
small. Using Drude model of free electron collisions {ly^ 
and ly^), permittivity and permeability read respectively 
e = 1 - + ii^e)]~^ and ^ = 1 - lu'^[lu{uj + 

where ujp and ujm are electric and magnetic plasma fre- 
quencies IS:]. For linearly polarized EM pulses (E = 
xS, H — yH) propagating along the z axis (k||z) Eqs. 
(1-2) become scalar ones. Applying paraxial approxi- 
mation (fc >> V) and neglecting higher order z deriva- 
tives Eqs. (1-2) can be reduced to generalized Ginzburg- 
Landau equation in NIMs (GLENIM) for electric field 



;dE 



K' d^E I 1 A , F -I- j'tttMEHI 



2fe 



E 



(3) 



and the equivalent magnetic field GLENIM obtained us- 
ing medium impedance Z = (/^/e)^/-^ — E/H. The fol- 
lowing renormalisation is used: ujpT t, {ujp/c)(^ — ^ z, 
uj/ojp vj, and (c/wp)^A^ —> = d'^E/dx^ + 

d^E/dy^. The rcnormalized frequency w runs in Fig.l 
till ujm/'^p — 0.8 |J^]. In this range the group velocity Vg 
is always positive and the phase velocity is negative. The 
group velocity dispersion (GVD) /c" = —{nvj)~^'^\W] is 
anomalous, thus, negative till vj — 0.7, changing sign af- 
ter (see also Fig. 2). As many materials, NIMs also can 
exhibit cubic Kerr nonlinearity which may be saturated 
, [l3l ■ In order to prevent pulse collapse the cubic non- 
linearity is usually saturated by a quintic one with the op- 
positesigne„i = [e^^) +ief^)\E\''-{ef) +isf^)\E\'^ 

and = (m1'^ + ^^if'^)\H\■' - (^f'^ + ^^if'^)\H\\ 

The sign of the real part of these indices determines 
the focusing/defocusing properties of the medium while 
imaginary parts are related to the nonlinear loss and gain. 
Nonlinear permittivity and permeability in synergy de- 
fine the cubic susceptibility x^^'' = s^^^ Z -\- ^'^^'^ /Z^ as well 
as the quintic one x^^-* = e''^^ Z + ^^^'^ / Z^ . Corresponding 
relations for magnetic susceptibility xh can be obtained 
using the impedance Z . Therefore, the system of coupled 
Eqs. (1-2) is reduced to a single GLENIM 



oz 



i{5E + 13 



a'E 



aj_Aj_E + K \E\^ E + v \Ef E = 
- eil^E + 'q\E\^ E + \Ef E) 



(4) 



with the renormalisation (2/ \k"\Y/'^t t, A±/2 \ n\ VJ 
Aj^, and ( 



wj^Yl'^E E. In reality, every ma- 
terial exhibits losses. The dissipative parameter of lin- 
ear loss 8 = —'Uj'^\e\ji\l1n is always negative insuring 



background stability since in NIMs %e[i\ < [2. The 
spectral filtering term is positive taking into account 
that /3 = '^\yV\l(nvj\k"\) > since 5[I^] < 0, as 
well as the energy diffusion term with 6 = 2ki/k^ > 

(3) (3) 

where ki > 0. Positive f] = Xi / Xr corresponds to 
the cubic gain compensating linear and nonlinear losses. 



thus, 
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^(3) 



must be negative [llj. The pa- 
rameter (T = sgn{—k") is either positive for anomalous 
GVD (AGVD) or negative for normal GVD (NGVD) 
(see Fig.l). In NIMs the diffraction has a negative sign 
since the coefficient a±_ = sgn{n) < 0. Therefore, bright 
spatial solitons exist only if the cubic nonlinearity is of 
the same sign as diffraction, thus, k — sgn{xr ) < 0. 
Consequently, the quintic nonlinearity has to be positive 



(3) 
Xr 



> 0. 



In order to generate bright spa- 

tiotemporal solitons so-called light bullets both diffrac- 
tion and dispersion must be compensated by saturating 
nonlinearity [ll[. Hence, both linear and nonlinear ef- 
fects need to have the same sign. As a consequence, light 
bullets can be generated only in NIMs with NGVD in the 
range 0.7 < va < 0.8 for the choice of dissipative param- 
eters P and 6 as in Fig. 2. The lack of space limits our 
studies here to temporal solitons described by temporal 
GLENIM (TGLENIM) corresponding to Eq. (4) without 
spatial second derivatives. The generation, propagation, 
and stability of light bullets in NIMs are investigated 
elsewhere. The obtained nonintegrable TGLENIM can 
be solved only numerically. However, some analytical 
approach even thought approximate is needed in order 
to have a better physical inside. The approximate reso- 
lution of TGLENIM is done using variation method ex- 
tended to dissipative systems fll'] . Although it is unable 
to account structural change of pulse profile, this method 
can serve as a guideline for simulation. A trial function 
E = AA^ exp [t^{iC - 0.5T-2)T-2 + it/;] is chosen with 
amplitude A{z), pulse temporal width T{z), wave front 
curvature C{z), and phase V'(^)- Using renormalisation 
of the propagation variable z/T^ — > z and linear loss 
parameter So = \S\ « 1, with = 2^/^ (4/3)^/^, 
and = (3/2)'^^^ 2^^^'^, these functions are optimized 
in order to yield corresponding Euler-Lagrange ordinary 
differential equations 



dA 

dz 



{So 



^A^ + 27^4 



T2 



2aC)A, (5) 



^ = {4aC - 11 A^ - 7A4)r + I - , (6) 

dz T 



dz 



a 



A^ 



A^ 



''rp2 '^^2 



4Ac 



dip 0" 5 ,9 

— = ^ + -kA'^ 

dz T2 2 



2iyA^ + 2/3C. 



(7) 



(8) 



Solving Eqs. (5-7) with zero z derivatives, a double 
steady-state solution with amplitudes A^ and A^ is 
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obtained Al ^ i0>^^-^r,)±^i^^^-^fW,^P^^) _ ^he 

± 2(/3i'CT— 37) 

family of solitons in conservative systems, reduces into 
either double or simple solution for each set of dissi- 
pative parameters. A double solution (A_ > A^) ex- 
ists in the (77, 7) domain between the parabola (/Jkct — 
+ 8(5o(37 - I3va) = and straight hne A_ 1 
for NGVD and AGVD in Fig.3. Above A_ = 1 per- 
sists only A+. The ratio between dissipative parame- 
ters /3 and 5 (given in Fig. 2) is kept same for NGVD 
and AGVD, in order to stress t he s imilarity of do- 
mains. The beam power P = i^/^TlA^T is no more 
conserved in dissipative systems However, the 

temporal width T = A~^{aK, + avA^)~'^/'^ and the 
power depend, up to v = max{|/3| , , |?7| , I7I}, only 
on the amplitude. The striking difference from the con- 
servative systems is the nonzero wave front curvature 
C = -Al {(3k -r](j + {I3v - -ia)Al) 4"!. To be soliton 
a steady state solution has to be stable. Our stability 
criterion based on the method of Lyapunov's exponents 
has to be extended to NIMs in order to check the stability 
of each steady-state. The solution of Eqs. (5-7) is stable 
if and only if the real parts of the solutions A of the 
equation (A"^ + aiA^ + a2A + as) = are all nonpositive 
[ll|. The stability criterion for TGLENIM is satisfied 
when coefficients a fulfill Hurwitz conditions 

a2 = 4kA'^(k -t- vA^) > 0, 
as = 2A^{k + uA^f[2A^{l3v(T - 87) -I- (/Jkct - 47/)] > 0, 
04 = 4(k + vA^)[A^{l2-iHi - Srjiy) + A^{13t]k ~ 8Soi^)] > 0, 

(9) 

where = aia2 — as taking into account that ai = 
A'^{APiya - 67) -I- A'^{(3k(t - A-q). The steady-state solu- 
tion A- satisfies this stability criterion only between the 
parabola and curve s corresponding to a4 = in Fig.3. 
The solution A^ is everywhere unstable. However, the 
stability of the solution A- is only a prerequisite to ob- 
tain a soliton after a self-organizing evolution. Indeed, 
an input pulse chosen in the established stable domain 
corresponds to a point on the upper stable branch of an- 
alytically obtained bifurcation curve v for power P ver- 
sus control parameter 7 in Fig. 4. The curve v is only 
a good approximation of the numerically obtained exact 
bifurcation curve n composed of solitonic attractors for 
different parameters. The lower unstable branch corre- 
sponds to unstable solutions A+ . A stable steady state on 
the curve v taken as the input for numerical simulations 
evolves toward the soliton on the curve n. The domain 
of stability is checked point by point in order to confirm 
by numerical simulations that corresponding inputs al- 
ways lead to a soliton. Indeed, the numerically obtained 
domain of stability limited by the curve s„ in Fig.3 is 
even slightly larger than the analytical one. Numerical 



simulations confirm that a soliton is propagating with, 
analytically predicted, nonzero wave front curvature. As 
a consequence, the dispersion is overcompensated by sat- 
urating nonlinearity leading to the collapse [ll|. How- 
ever, the collapse is prevented by losses equilibrated in 
turn by gain. Therefore, the self-organization of dissipa- 
tionless solitons is based on crosscompensation. During 
their propagation such solitons render the medium effec- 
tively dissipationless. In order to check soliton robustness 
its amplitude is increased 40% at each z systematically 
3000 times (see Fig. 5). Such tremendous perturbations 
drastically increase the pulse temporal width and am- 
plitude. However, far from being annihilated the pulse, 
remaining in attraction domain without losses and dis- 
persion, maintains its new form. Indeed, to each increase 
the self-organized system reacts lowering its amplitude 
prj . Perturbations arrested the soliton recovers its ini- 
tial shape after about 400 steps in z giving evidence of 
astonishing robustness. Perturbations may correspond to 
monochromatic EM pulses injected in such novel active 
medium composed of dissipationless solitons and NIM 
(SOLINIM). Therefore, it seems that SOLINIM behaves 
for other pulses as effectively dissipationless. 

In conclusion, SOLINIM may be considered as a novel 
very promising active composite medium due to the syn- 
ergy between dissipationless solitons and NIM. In or- 
der to describe such SOLINIM systems, newly estab- 
lished coupled nonparaxial Helmholtz equations for elec- 
tric and magnetic fields are reduced to novel paraxial 
complex cubic-quintic Ginzburg-Landau equation. In 
contrast with ordinary media, bright dissipationless light 
bullets may propagate only in NGVD NIMs. However, 
bright dissipationless temporal solitons can be generated 
in NIMs with both NGVD and AGVD. A stabihty cri- 
terion is established. Choosing as input a steady-state 
with dissipative parameters from obtained stability do- 
mains for either NGVD or AGVD NIMs, a self-organized 
propagation always results in generation of extremely ro- 
bust dissipationless soliton. Such solitons during their 
propagation render NIM effectively dissipationless. We 
hope that in SOLINIMs the practical realization of pe- 
culiar effects like the cloaking and the superresolution 
will be no more prevented. 
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Figure captions 

Fig. 1. Group velocity Vg, GVD, and phase velocity Vph- 
Fig. 2. Dissipative parameters /3 and 5 versus frequency 
zu. 

Fig. 3. AGVD and NGVD domains of stable solutions 
A-. 

Fig. 4. Upper stable and lower unstable branches of vari- 
ational V curve and numerical n curve. 
Fig. 5. Dissipatiorilcss soliton resisting 40% increase of 
the amplitude for each z from z = 1000 to 2; = 4000. 
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